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Abstract 


ote- 


We  consider  &  class  of  two-point  boundary  val^e  problems^.  _ 
systems  of  the  form(>5)  =  f(x,y,t,e),  ey'rfi^^(x,t,e)  + /g^j(x,t,e)y  where 


or 


the  matrix^j^yx, t , 0)  has  a  hyperbolic  splitting  with^a" fixed  number  of 

stable  and  unstable  eigenvalues.  Solutions  to  such  boundary  value 
problems  can  then  be  expected  to  have  boundary  layer  behavior  near 

both  endpoints  in  the  limIt~^y<^-4i-^jrT7eobtain  uniform  asymptotic  re¬ 
presentations  of  solutions .  ^Our 'analysis  shows,  in  particular,  how  to 


determine  the  reduced  order  boundary  value  problem  satisfied  by  the 
limiting  interior  solution.  Orthogonal  matrix  methods  are  used  to 
determine  this  reduced  problem  and  appropriate  boundary  layer 
corrections  in  a  computationally  effective  manner.  Numerical  experi¬ 
ments  with  model  problems  illustrate  the  possibility  of  multiple  solu' 
tions  and  show  how  our  asymptotic  results  can  be  used  in  combination 
with  the  COLSYS  code  for  solving  two-point  problems  via  collocation. 


Initial  value  problems  for  stiff  systems  of  ordinary  differential 
equations  are  now  considered  to  be  relatively  tractable  numerically 
(cf.  Enright  ct  al.  (1975)),  though  still  the  subject  of  current  re¬ 
search  and  difficulty  (cf.,  e.g.,  Miranker  (1979)  and  Kreiss  (1978)). 
Codes  appropriate  for  two-point  problems  whose  solutions  involve  rapid 
changes  near  the  endpoints  are  not  readily  available,  though  these 
problems  occur  in  many  applications.  Wo  must  study  them  as  a  prelude 
to  analyzing  more  difficult  problems  with  interior  layers  of  rapid 
transition  at  unknown  locations. 


*This  work  was  supported  in  part  by  the  Air  Force  Office  of 
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Specifically,  let  us  consider  singular  perturbation  problems  of 
the  form 

(1)  x  =  f(x,y,t,e)  , 

(2)  Gy  =  g(x,y,t,C)  =  g^x.t.C)  +  g2(x,t,e)y  , 

(3)  A(x(0) ,y (0) ,C)  =  a  (x(0),e)  +  a2 (x(0) ,c)y(0)  =  0  , 

(4)  B(x(l),y(l),c)  =  b  (x(l),e)  +  b2 <x(l) ,c) y (1)  =  0  , 

where  x  and  y  are  vectors  of  dimensions  m  and  n,  respectively,  and 
where  the  n><n  matrix  g2(x,t,0)  has  k  >  0  (strictly)  stable  and  n-k  >  0 

(strictly)  unstable  eigenvalues  for  all  x  and  for  0  <_  t  <_  1.  We  shall 
assume  that  there  are  q  >  k  initial  conditions  (3)  and  r  =  m  +  n  -  q 
>  n-k  terminal  conditions  (4),  and  that  the  coefficients  in  (l)-(4) 
have  asymptotic  expansions  in  e  with  all  terms  being  smooth  functions 
of  their  remaining  arguments.  Limiting  solutions  of  ( 1 ) — (4)  will  be 
sought  as  the  small  positive  parameter  C  tends  toward  zero. 

O'Malley  (1980)  discusses  the  behavior  of  solutions  to  more  non¬ 
linear  problems.  We've  simplified  that  treatment  by  supposing  that 
the  (potentially)  stiff  differential  system  (2)  and  the  separated 
boundary  conditions  (3) -(4)  are  linear  in  the  "fast"  variable  y. 
Vasil'eva  and  Butuzoy  (1978)  and  O'Malley  and  Flaherty  (1980)  discuss 
problems  where  g2(x,t,0)  has  a  nullspace,  while  Miranker  and  Wahba 

(1976)  and  Kreiss  (1979)  discuss  problems  where  g  (x,t,0)  has  purely 
imaginary  eigenvalues.  The  numerical  treatment  of  certain  interior 
and  shock  layer  problems  is  considered  in  Kreiss  (197S) ,  Hemker  (1977) 
Osher  (1980),  and  elsewhere. 

For  the  special  problem  (l)-(4),  we  can  expect  bounded  solutions 
to  have  the  limiting  form 

x(t,G)  =  XQ(t)  +  0(G)  , 

(5) 

y  (t,c)  =  Y  ( t)  +  |J  (t/G)  +  V  ((l-t)/G)  +  0(G)  , 

0  0  o 

although  unbounded  solutions  to  such  problems  are  also  possible  (cf. 
O'Malley  (1970)  and  Ferguson  (1975)).  Indeed,  impulsive  endpoint  be¬ 
havior  is  important  in  applications  to  control  and  systems  theory  (cf. 
O'Malley  (1978)).  Here,  we'll  ask  that  Mq(t)  and  vq(cj)  decay  to  zero 

as  T  and  0,  respectively,  tend  to  infinity.  Then,  p  will  provide  a 

rapid  transition  from  y(0,0)  to  YQ(0)  near  t  =  0,  and  VQ  will  allow 

analogous  nonuniform  convergence  as  G  -►  0  at  the  terminal  endpoint. 

The  limiting  solution,  (Xq,Yq),  within  (0,1)  will  necessarily  satisfy 

the  limiting  system  XQ  =  f (Xn» Y^, t,0) ,  g(X^,Y0,t,0)  =  0,  so 

(6)  YQ(t)  =  G(XQ,t)  =  -g"1(X0,t,0)g1(X0,t,0)  ,  !  TI  I  i 
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and  XQ  must  satisfy  the  mth  order  nonlinear  (and  non-stiff)  system 

(7)  XQ  =  F(XQ,t)  =  f (X0,G(X0,tO  ,t,0)  . 

We  need  to  specify  an  appropriate  set  of  m  boundary  conditions 
for  our  limiting  solution  Xg(t)  =  x(t,0)  within  (0,1).  Thus,  we  seek 

a  "cancellation  law"  (cf.  Wasow  (1944)  and  Harris  (1973))  which 
selects  a  combination  of  q-k  of  the  limiting  initial  conditions  (3) 
and  of  r-n+k  of  the  terminal  conditions  (4)  to  be  satisfied  by  XQ. 

The  limiting  solution  (Xq»Yq)  will  not  generally  satisfy  the  m  limiting 

boundary  conditions,  i.e.  AQ(X0(0))  =  A(XQ(0) ,YQ(0) ,0)  and  B^(XQ(1))  = 

B(Xq(1) ,Yq(1) ,0)  will  not  vanish.  We  will,  however,  find  m  boundary 

conditions  for  XQ  of  the  form 

(8)  <J>(X0(0))  =  Z0+(V°))A0(X0<0>>  =  0  ' 

(9)  f(X0(l))  =  Z1_<X0(1))B1(X0(1) )  =  0  , 

for  appropriate  matrices  ZQ+  and  Z^_  with  ranks  q-k  and  r-n+k,  respec¬ 
tively.  The  nonlinear  "reduced"  two-point  problem  (7) -(9)  can  be 
expected  to  have  no  solution,  a  unique  solution,  or  many  solutions. 

To  begin  our  humerical  discussion,  recall  that  the  Schur  decompo¬ 
sition  (cf.,  e.g.,  Stewart  (1973))  guarantees  the  existence  of  an 
orthogonal  matrix  E  such  that 

T  ( x , t )  U(x,t) 

(10)  q?(x,t,0)E(x,t)  =  E(x,t)  "  Q  T 

where  T_  is  k*k  and  upper  triangular  with  the  stable  eigenvalues  of  q^ 

and  where  T+  is  upper  triangular  with  n-k  unstable  eigenvalues.  The 

matrix  E  can  be  determined  numerically  through  the  QR  algorithm,  i.e. 
through  a  sequence  of  Householder  transformations.  We  actually  need 
much  less  than  this,  and  would  prefer  a  procedure  which  only  block 
triangularized  q^.  Partitioning 

(11)  E  =  [E_  E+]  , 

after  its  kth  column,  note  that  E  (E+)  spans  the  stable  (unstable) 
eigenspace  of  q^,  so  P  =  E_E*  and  Q  =  E+E^  are  complementary  projec¬ 
tions  onto  these  eigenspaces. 

In  order  to  have  a  limiting  solution  of  the  form  (4),  the  initial 
layer  correction  term  Hq(t)  must  be  a  decaying  solution  of 


(13) 


U0<T)  =  exp(G0T)viQ(0) 


i.e. 


will  decay  provided  UQ(0)  lies  in  the  stable  eigenspace  of  GQ, 

Mo(0)  "  Eo-co  where  V  =  E-{X0(0)'0)* 

Since  we  are  seeking  solutions  of  the  form  (5) ,  our  representa¬ 
tions  (6)  and  (13)  imply  that  the  q  limiting  initial  conditions  (3) 

take  the  form  A(x(0,0) ,y(0,0) ,0)  =  A(X  (0),Y  (0)  +  E  c  ,0)  =  A  (X  (0)) 

O  .j.  \j  0“  0  o  u 

+  a2(Xg(0)  '0)E0_cq  =  0,  so  cQ  =  ~(a2oEo  ^  Ao  where  tlle  da99er  repre¬ 
sents  a  matrix  pseudoinverse.  Assuming  that  a20Eo  ^as  max^ma^ 
rank  k,  there  will  be  an  orthogonal  matrix  ZQ  =  (Z^  Z^+)  *  which  re¬ 

duces  it  to  its  row  echelon  form,  i.e.  Z  a  E  =  (V’  0)’  where  V  is 

0  20  0—  0  0 

k*k  and  nonsingular.  After  multiplying  by  ZQ,  the  first  k  rows  of  the 

product  provide  the  unknown  (but,  due  to  the  flexibility  in  selecting 
ZQ,  not  necessarily  unique)  vector 

1141  “o'01  ■  -Eo-vo‘VAo  • 

needed  to  specify  HQ(t).  The  remaining  q-k  equations  of  the  product 

provide  the  initial  conditions  (8)  needed  to  define  the  reduced  prob¬ 
lem  for  XQ (t) . 

In  analogous  fashion,  if  we  suppose  that  b^E  =  b2(XQ(l),0)  • 
E+(Xq(1),1)  has  its  maximal  rank  n-k,  and  let  Z^  =  (Z|+  Zj  )  ’  be  such 
that  z1b2lEl+  =  *Vi  '  *s  r°WTeduced,  the  decaying  terminal  boundary 
layer  correction  term 

(15)  VQ(o)  .  -exp(g2(Xo(l) ,l,0)o)E1+v“1Z1+Bo  , 

becomes  (usually  nonuniquely)  specified  as  does  the  set  of  terminal 
conditions  (3)  for  X^ft). 

Because  the  orthogonal  matrices  EQ,  ZQ,  ,  and  Z^  depend  on  the 

corresponding  m  vectors  Xq(0)  =  x(0,0)  and  X^ll)  =  x(l,0),  we  shall 

usually  have  to  use  Newton's  method  to  actually  determine  them  when¬ 
ever  g  (x,t,0),  a2(x,0),  or  b2(x,0)  actually  depend  on  the  x  variable 

at  t  =  0  or  1.  [Alternatively,  unless  the  boundary  conditions  (3)  and 
(4)  directly  prescribe  or  determine  x(0,e)  or  x(l,c),  a  shooting 
technique  (cf.  Keller  (1968)),  which  involves  guessing  the  endvalue 
X^ (0)  or  Xq(D  of  the  slow  x(t,0)  vector  and  integrating  XQ  to  the 

opposite  endpoint,  could  be  attempted.  Such  a  procedure  would  be  far 
better  than  shooting  for  the  full  vector  (x,y),  an  approach  which  does 
not  work  well  for  problems  with  boundary  layer  behavior  at  both 


endpoints.]  Newton's  method  reouires  a  satisfactory  initial  guess 
0 

X  (t)  for  X  (t) .  We  would  proceed  to  determine  successive  iterates 
V  0 

X  (t)  until  we  achieved  numerical  convergence.  At  the  vth  step,  we 

\> 

would  calculate  an  approximation  E^  to  E(XQ(p),p,0)  for  p  =  0  and  1  by 
applying  the  QR  procedure  to  g2 (XV(p) ,p,0) .  One  QR  step  with  the 

matrix  ZV  will  reduce  each  of  a_(XV(0),0)E  (XV(0),0)  and  b„(XV(l),0)  • 

V  P  2  V2  v+1 

E+(X  (1),1)  to  row  echelon  form.  A  linear  problem  for  6x  (t)  =  X  (t) 

v 

-  X  (t)  is  thereby  determined,  presuming  we  can  neglect  the  influence 
of  derivatives  of  ZQ  and  Z  on  the  selection  of  the  boundary 
conditions.  ‘J 


To  illustrate  such  problems  numerically,  we  have  experimented 
with  the  simple  harmonic  oscillator  system 

(16)  x  =  1  -  x  ,  cy^^  =  y2  ,  ey2  =  a2(x)y1  +  6(x)  , 

made  nonlinear  since  the  coefficients  a(x)  =  1  +  2x  and  B(x)  =  8x(l-x) 
depend  on  the  slow  vector  x.  The  three  linear  boundary  conditions  are 


(17)  x(0)  +  yx (0)  =  0  ,  -bx (0)  +  y2(0)  =  0  ,  and  x(l)  +  y.(l)  =  0. 

Since  g2  has  one  positive  and  one  negative  eigenvalue  whenever  a  is 

nonzero,  we  can  expect  the  limiting  solution  within  (0,1)  to  be  deter¬ 
mined  by  a  cancellation  law  which  retains  only  one  initial  condition. 
Thus,  the  reduced  problem  will  consist  of  the  limiting  system 


(18) 


xo  ■  1  -  *0  ■ 


\o  ■  0  •  “  <VV10  +  BIV  ■  0  • 


and  some  combination  of  X„(0)  +  Y,  (0)  and  -bx  (0)  set  to  zero. 

0  10  0 

The  E  matrix  for  this  example  is  explicitly  determined  as 
^  (  i  .  .  n 


E(x,t)  = 


+  a2  (x) 


1  |a(x) 

(- |ot  (x)  I  1 


,  and  we  can  take  ZQ  to  be  any 


multiple  of  E' (XQ(0) ,0) .  Thus,  the  initial  condition  (8)  appropriate 
for  the  reduced  system  (19)  is 


(20)  |ct(X0(0) )  |  (X  (0)  +  Y1Q(0))  -  bXQ(0)  =  0  . 

Rewriting  this  as  a  function  of  xfi  =  XQ(0)  alone,  xfl  must  satisfy  a 

cubic  equation  with  the  three  roots  x  =  0  and  ^]b  sgn  a  -  6 
2  1/2  0  40 

+  ((b  sgn  -  4)  +  48)  ]  selected  so  that  sgn  =  +1  according  to 

the  sign  of  1  +  2xQ.  (For  b  =  2,  the  corresponding  initial  values  are 

xQ  =  0,  0.803,  and  -4.29.)  Values  of  xR  must  be  avoided  for  which  the 

corresponding  a(Xg)  =  3  +  2(xQ  -  l)e  has  a  zero  within  0  <  t  <  1. 


5 


Since  XQ(0)  +  V^g(O)  =  -bX  (0)/(l  ■■  2XQ(0))  neither  of  the  separate 

limiting  initial  conditions  h^lX^U'))  =  0  will  be  satisfied  unless 

bXg(O)  =  0.  Only  then,  will  the  initial  boundary  layer  term  \i Q (t )  be 

trivial.  Unless  Y^fl)  =  -X  (1)  =  -1  _  (x^  -  l)e  \  the  terminal 

boundary  layer  term  VQ(a)  is  also  nontrivial. 

We  are  guaranteed  that  the  three  “outer”  solutions  <x0«^q)  s9 

obtained  provide  limiting  solutions  within  (0,1)  for  corresponding 
actual  solutions  in  the  limit  £  ■*■  0  (cf.  Hoppensteadt  (1971)).  More¬ 
over,  adding  the  corresponding  boundary  layer  corrections  iiQ(t/e)  and 

vQ( (l-t)/e)  (cf.  (13)  and  (16))  will  improve  the  outer  approximation 

YQ(t)  in  the  endpoint  boundary  layers. 

One  would  expect  our  approximation  (5)  to  be  a  good  first  guess 
to  provide  a  two-point  boundary  value  solver.  This  is  not  necessarily 
true,  however,  because  of  stiffness  as  c  +  0.  Our  experience,  instead, 
suggests  usinq  our  (then  crude)  asymptotic  solution  as  a  first  guess 
for  problems  where  £  is  only  moderately  small  and  then  gradually 
reducing  £.  This  scheme  allows  mesh  points  to  be  gradually  redistri¬ 
buted  in  boundary  layer  regions  of  rapid  transiton,  and  provides  good 
numerical  results. 

We  used  the  COLSYS  code  of  Ascher  et  al_.  (1979),  but  note  that 
more  specialized  methods  for  singular  perturbation  problems  might  also 
be  tried  (cf.  Flaherty  and  Mathon  (1980)).  We  calculated  three  dis¬ 
tinct  solutions  to  our  boundary  value  problem  for  b  =  0,  1,  2,  and  -2 
following  the  procedure:  The  asymptotic  solution  and  a  uniform  mesh 

-1  -1  -2 

was  used  as  an  initial  guess  for  E  =  10  .  The  E-sequence  10  ,10  , 

-4  -6 

10  ,  and  10  was  used,  wj  th  the  previously  computed  solution  and 

mesh  used  as  the  initial  approximation.  If  the  solution  failed  to 

converge  or  used  more  than  222  sub  intervals,  C  =  10  ^  and/or  C  =  10  ^ 
values  were  inserted  into  the  sequence.  Calculations  used  both  cubic 

*6 

and  quintic  splines,  and  error  tolerances  of  10  for  x  and  y  and 

—  3  —6  A 

either  10  or  10  for  y^.  In  general,  good  accuracy  was  obtained 

with  both  cubic  and  quintic  splines,  though  calculation  time  was 
generally  less  with  the  lower  order  polynomials.  Alternative  solu¬ 
tions  with  slightly  different  initial  quesses  took  much  more  computing 
time.  For  b  =  -2  and  XQ(0)  *  -2.80,  we  have  a  turning  point  since 

tt(X  (0.930))  =  0.  Actual  solutions  become  unbounded  there  as  E  -*■  0, 

0  -2 
though  for  moderate  values  of  E  (like  10  )  COLSYS  could  provide  a 

corresponding  solution  with  difficulty.  For  b  =  0,  and  Xg(0)  =  “3.5, 

a  has  a  zero  above  t  =  1.  This  forces  the  boundary  layer  jump  |y^0(0) 

-  y(l,0)|  to  be  large  (89.8  compared  to  0.4  and  0.6  for  the  other 
roots  XQ(0)  =  0  and  0.5),  yet  we  get  good  results  for  e  =  10  As  a 
typical  example,  we  display  our  results  (Figures  1,  2,  and  3)  for 


b  =  2  and  tabulate  our  results  (Table  1)  for  cubic  splines.  The 
column  headings  for  Table  1  ares  (1)  e,  (2)  the  number  of  subinter¬ 
vals  used  by  COLSYS  to  solve  the  problem  with  the  given  tolerances, 

(3)  and  (4)  the  number  of  subintervals  in  the  right  (left)  boundary 
layer,  defined  as  t  <  10c  ( t>  1  -  loe),  (5)  the  percentage  of  subin¬ 
tervals  in  the  boundary  layers,  (6)  time  to  perform  the  calculation, 
excluding  input/output,  and  (7)  accumulated  CP  time  on  an  IBM  3033. 
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